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Formation of sharp eccentric rings in debris disks 
with gas but without planets 

W. Lyra 1 ’ 2 ’ 3 & M. Kuchner 4 


‘Debris disks’ around young stars (analogues of the Kuiper Belt in 
our Solar System) show a variety of non-trivial structures attri- 
buted to planetary perturbations and used to constrain the pro- 
perties of those planets 1-3 . However, these analyses have largely 
ignored the fact that some debris disks are found to contain small 
quantities of gas 4-9 , a component that all such disks should contain 
at some level 10,11 . Several debris disks have been measured with a 
dust-to-gas ratio of about unity 4-9 , at which the effect of hydro- 
dynamics on the structure of the disk cannot be ignored 12,13 . Here 
we report linear and nonlinear modelling that shows that dust-gas 
interactions can produce some of the key patterns attributed to 
planets. We find a robust clumping instability that organizes the 
dust into narrow, eccentric rings, similar to the Fomalhaut debris 
disk 14 . The conclusion that such disks might contain planets is not 
necessarily required to explain these systems. 

Disks around young stars seem to pass through an evolutionary phase 
when the disk is optically thin and the dust-to-gas ratio s ranges from 0. 1 
to 10. The nearby stars (3 Pictoris 5,6,15-17 , HD32297 (ref. 7), 49 Ceti (ref. 4) 
and HD 21997 (ref. 9) all host dust disks resembling ordinary debris 
disks and also have stable circumstellar gas detected in molecular CO, 
Na i or other metal lines; the inferred mass of gas ranges from lunar 
masses to a few Earth masses (Supplementary Information). The gas in 
these disks is thought to be produced by planetesimals or dust grains 
themselves, by means of sublimation, photodesorption 10 or collisions 11 , 
processes that should occur in every debris disk at some level. 

Structures may form in these disks by a recently proposed instability 12,13 . 
Gas drag causes dust in a disk to concentrate at pressure maxima 18 ; how- 
ever, when the disk is optically thin to starlight, the gas is most probably 
primarily heated by the dust, by photoelectric heating. In this circum- 
stance, a concentration of dust that heats the gas creates a local pressure 
maximum that in turn can cause the dust to concentrate more. The result 
of this photoelectric instability could be that the dust clumps into rings or 
spiral patterns or other structures that could be detected by coronographic 
imaging or other methods. 

Indeed, images of debris disks and transitional disks show a range of 
asymmetries and other structures that call for explanation. Traditionally, 
explanations for these structures rely on planetary perturbers — a tantali- 
zing possibility. However, so far it has been difficult to prove that these 
patterns are clearly associated with exoplanets 19,20 . 

Previous investigations of hydrodynamical instabilities in debris disks 
neglected a crucial aspect of the dynamics: the momentum equations for 
the dust and gas. Equilibrium terminal velocities are assumed between 
time steps in the numerical solution, and the dust distribution is updated 
accordingly. The continuity equation for the gas is not solved; that is, the 
gas distribution is assumed to be time-independent, despite heating, 
cooling, and drag forces. Moreover, previous investigations considered 
only one-dimensional models, which can only investigate azimuthally 
symmetrical ring-like patterns. This limitation also left open the possibi- 
lity that, in higher dimensions, the power in the instability might collect in 
higher azimuthal wavenumbers, generating only unobservable clumps. 


We present simulations of the fully compressible problem, solving 
for the continuity, Navier- Stokes and energy equations for the gas, and 
the momentum equation for the dust. Gas and dust interact dynami- 
cally through a drag force, and thermally through photoelectric heating. 
These are parametrized by a dynamical coupling time r f and a thermal 
coupling time t t (Supplementary Information). The simulations are 
performed with the Pencil Code 21-24 , which solves the hydrodynamics 
on a grid. Two numerical models are presented: a three-dimensional 
box embedded in the disk that co- rotates with the flow at a fixed dis- 
tance from the star; and a two-dimensional global model of the disk in 
the inertial frame. In the former the dust is treated as a fluid, with a 
separate continuity equation. In the latter the dust is represented by 
discrete particles with position and velocities that are independent of 
the grid. 

We perform a stability analysis of the linearized system of equations 
that should help interpret the results of the simulations (Supplemen- 
tary Information). We plot in Fig. la-c the three solutions that show 
linear growth, as functions of e and n = kH , where k is the radial 
wavenumber and H is the gas scale height ( H = c s / s/yQ K > where c s 
is the sound speed, f2 K the Keplerian rotation frequency and y the 
adiabatic index). The friction time T f is assumed to be equal to l/Q K . 
The left and middle panels show the growth and damping rates. The 
right panels show the oscillation frequencies. There is no linear insta- 
bility for £ > 1 or n < 1. At low dust load and high wavenumber the 
three growing modes appear. The growing modes shown in Fig. la 
have zero oscillation frequency, characterizing a true instability. The 
two other growing solutions (Fig. lb, c) are overstabilities, given the 
associated non-zero oscillation frequencies. The pattern of larger 
growth rates at large n and low £ invites us to take £ = sn 2 as character- 
istic variable and to explore the behaviour of 1 . The solutions in this 
approximation are plotted in Fig. If, g^ The instability (red) has a 
growth rate of roughly 0.26f2 K for all £. The overstability (yellow) 
reaches an asymptotic growth rate of Q K / 2, at ever-growing oscillation 
frequencies. Damped oscillations (blue) occur at a frequency close to 
the epicyclic frequency. 

Whereas the inviscid solution has growth even for very small wave- 
lengths, viscosity will cap power at this regime, leading to a finite fastest- 
growing mode (Supplementary Information), which we reproduce 
numerically (Fig. lh). Although there is no linear growth for s > 1, we 
show that there exists nonlinear growth for s = 1. We show in Fig. li the 
time evolution of the maximum dust surface density Z d (normalized by 
its initial value, 27 0 ). A qualitative change in the behaviour of the system 
(a bifurcation) occurs when the noise amplitude of the initial velocity 
(Wrms) is raised far enough, as expected from nonlinear instabilities 25,26 . 
We emphasize this result because, depending on the abundance of H 2 , 
the range of s in debris disks spans both the linear and nonlinear regimes. 
The parameter space of t t and Tf is explored in one-dimensional models 
in Supplementary Information, showing robustness. 

In Fig. 2 we show the linear development and saturation of the 
photoelectric instability in a vertically stratified local box of size 
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Figure 1 | Linear analysis of the axisymmetric modes of the photoelectric 
instability. Solutions for axisymmetric perturbations i/f = \j/ exp(st + ikx), 
where \j/ is a small amplitude, x is the radial coordinate in the local Cartesian co- 
rotating frame, k is the radial wavenumber, t is time and s is the complex 
frequency. Positive real s means that a perturbation grows, negative s indicates 
that a perturbation is damped, and imaginary s represents oscillations. 
Solutions are for a = 0, Tf = l/f2 K and t t = 0. a-e, The five solutions as 
functions of n = kH and e. Solutions a-c show linear growth. Growth is 
restricted to the region with low dust-to-gas ratio (e < 1 ) and high wavenumber 
(n > 1). The growing modes in b and c have non-zero oscillation frequencies, 
characterizing an overstability. Conversely, solution a is a true instability, 
d, e, Solutions that correspond to damped oscillations through most of the 
parameter space. In a small region (high dust-to-gas ratio and high frequency), 


modes are exponentially damped without oscillating, f, g, Growth rate (f) and 
oscillation frequency (g). Using f = sn 2 and taking the limit f 1 permits better 
visualization of the three behaviours: true instability (red), overstability 
(yellow) and damped oscillations (blue). The other two solutions are the 
complex conjugate of the oscillating solutions and are not shown, h, Growth 
rates. When viscosity is considered (a = 10 -2 in this example), power is capped 
at high wavenumber, leading to a finite most-unstable wavelength. The figure 
shows the analytical prediction of the linear instability growth in this case 
(Supplementary Information) compared to the growth rates measured 
numerically. The overall agreement is excellent. The growth rates are only very 
slightly underestimated, i, Nonlinear growth. Although there is no linear 
instability for £ = 1, growth occurs when the amplitude of the initial 
perturbation ( u rms ) is increased, a hallmark of nonlinear instability. 


(1 X 1 X 0.6)H and resolution 255 X 256 X 128. The dust and gas are 
initialized in equilibrium (Supplementary Information). The dust-to- 
gas ratio is given by loge = —0.75, so that there is linear instability, and 
viscosity v = otc s H is applied as a = 10~ 4 (where a is a dimensionless 
parameter 27 ). The initial noise is u rnis /c s = 10 -2 . Figure 2a shows the 
dust density p^ in the x-z plane, and Fig. 2b that in the x-y plane, both 
at 100 orbits (the orbital period is T orb = 2n/Q K ). Figure 2c shows the 
one-dimensional x- dependent vertical and azimuthal average against 
time. Through photoelectric heating, pressure maxima are generated at 
the locations where dust concentrates, which in turn attract more dust 
by means of the drag force. There is no hint of unstable short- wavelength 


(less than H ) non- axisymmetric modes: the instability generates stripes. 
The simulation also shows that stratification does not quench the insta- 
bility. Figure 2d shows a plot of the maximum dust density against time, 
achieving saturation and steady state at about 70 orbits. 

We consider now a two-dimensional global model. The resulting 
flow, in the r-0 plane (r is radius and (j) is azimuth), is shown in Fig. 3a-c 
at selected snapshots. The flow develops into a dynamic system of nar- 
row rings. Whereas some of the rings break into arcs, some maintain 
axisymmetry for the whole timespan of the simulation. It is also observed 
that some arcs later re-form into rings. We check that, in the absence 
of the drag force back- reaction, the system does not develop rings 
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Figure 2 | Growth and saturation of the photoelectric instability. In this 
three-dimensional stratified local box with linearized Keplerian shear, the main 
source of heating is photoelectric. The equilibrium in the radial direction is 
between stellar gravity, Coriolis force and centrifugal force. In the vertical 
direction the equilibrium for the gas is hydrostatic, between stellar gravity, 
pressure and the drag- force back-reaction. To provide a stable stratification, an 
extra pressure pb = pc£ is added, where c h is a sound speed associated with a 
background temperature. For the dust, a steady state is established between 
gravity, diffusion and drag force. The dust continually falls to the midplane but 
is diffused upwards. The diffusion is applied only in z, mimicking turbulent 
diffusion that is in general anisotropic, a, x-z cut aty = 0 at 100 orbits. The 
instability concentrates dust in a preferred wavelength. The resulting structures 
have stable stratification, b, x-y cut at the midplane z = 0 at 100 orbits. No 


non-axisymmetric instability is observed, and the dust forms stripes, c, Time 
evolution of the vertically and azimuthally averaged density, showing the 
formation of well-defined rings, d, Time evolution of the maximum dust 
density. The instability saturates at about 70 orbits in this case. The slowdown 
compared with the growth rate Q K f 2 predicted in Fig. 1 is because of the use of 
viscosity, and the background pressure needed for the stratification. The 
dimensionless parameter f$ = y(c h /c s ) 2 measures the strength of this term, 
e, Maximum growth rate, showing that linear instability exists as long as < 1. 
The maximum growth rate decreases smoothly from Q K /2 for ft = 0, to zero for 
ft = 1. f, The structure formed in the dust density at t = 50 (about eight orbits) 
for different values of f$. At moderate values, growth still occurs at a significant 
fraction of the dynamical time. The run shown in a-d used p = 0.5. 


(Supplementary Information). We also check that when the conditions 
for the streaming instability 24 are considered, the photoelectric insta- 
bility dominates (Supplementary Information). 


A development of the model is that some of the rings start to oscil- 
late, seeming eccentric. These oscillations are epicycles in the orbital 
plane, with a period equalling the Keplerian, corresponding to the free 



Figure 3 | Sharp eccentric rings, a-c, Snapshots 
of the dust density in a two-dimensional global disk 
in polar coordinates, at 20 orbits (a), 40 orbits 
(b) and 60 orbits (c). The photoelectric instability 
initially concentrates the dust axisymmetrically 
into rings, at a preferred wavelength. As the 
simulation proceeds, some rings maintain the 
axisymmetry and others break into arcs. Some arcs 
rearrange into rings at later times, such as those at 
r = 0.6 and r = 1.0 between b and c. Although 
mostly axisymmetric, some rings seem to oscillate, 
appearing off-centred or eccentric, d, We measure 
the azimuthal spectral power of the density shown 
in c, as a function of radius. Modes from m = 0 to 
m = 3 are shown, where m is the azimuthal 
wavenumber, e, Although the ring at r = 1.5 has 
m = 0 as the more prominent mode, we show that a 
circle (black dotted line) is not a good fit. An ellipse 
of eccentricity e = 0.03 (red dotted line) is a better 
fit, although still falling short of accurately 
describing its shape. The black and red diamonds 
are the centre of the circle (the star) and the centre 
of the ellipse (a focal distance away from the star), 
respectively. 
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oscillations in the right-hand side of Fig. la-c. We check (Supplemen- 
tary Information) that they correspond to eigenvectors for which 
u = v; that is, gas and dust velocities coinciding. For this mode, the 
drag force and back- reaction are cancelled. So, for maintaining the 
eccentricity, this mode is being selected from among the other modes 
in the spectrum. This is naturally expected when the dust-to-gas ratio 
is very high. For 8^1, the gas is strongly coupled to the dust, cancel- 
ling the gas-dust drift velocity in the same way that if 1 does in the 
opposite way, by strongly coupling the dust to the gas. In this configu- 
ration, the freely oscillating epicyclic modes can be selected. 

We plot in Fig. 3e one of the oscillating rings, showing that its shape 
is better fitted by an ellipse (red dotted line) than by a circle (black 
dotted line). The eccentricity is 0.03, which is close to the eccentricity 
found 28 for the ring around HD 61005 ( e = 0.045 ± 0.015). We also 
notice that some of the clumps in Fig. 3 should become very bright in 
reflected light, because they have dust enhancements of an order of 
magnitude. In conclusion, the proposed photoelectric instability pro- 
vides simple and plausible explanations for rings in debris disks, their 
eccentricities, and bright moving sources in reflected light. 

Recent work 29 suggests that the ring around Fomalhaut is confined 
by a pair of shepherding terrestrial-mass planets, below the current 
detection limits. Detection of gas around the ring would be a way to 
distinguish that situation from the one we propose. At present, only 
upper limits on the amount of gas in the Fomalhaut system exist 30 ; 
however, they are relatively insensitive because they probe CO emis- 
sion, and CO could easily be dissociated around this early A-type star. 
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